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Abstract. Following the seminal work of F. Bouchut on zero pressure gas dynamics, extensively 
used for gas particle-flows, the present contribution investigates quadrature-based velocity moments 
models for kinetic equations in the framework of the infinite Knudsen number limit, that is, for dilute 
clouds of small particles where the collision or coalescence probability asymptotically approaches zero. 
Such models define a hierarchy based on the number of moments and associated quadrature nodes, 
the first level of which leads to pressureless gas dynamics. We focus in particular on the four moment 
model where the flux closure is provided by a two-node quadrature in the velocity phase space and 
provides the right framework for studying both smooth and singular solutions. The link with both 
the kinetic underlying equation as well as with zero pressure gas dynamics, i.e. the dynamics at the 
frontier of the moment space of order four, is provided. We define the notion of measure solutions 
and characterize the mathematical structure of the resulting system of four PDEs. We exhibit a 
family of entropies and entropy fluxes and define the notion of entropic solution. We study the 
Riemann problem and provide entropic solutions in particular cases. This leads to a rigorous link 
with the possibility of the system of macroscopic PDEs to allow particle trajectory crossing (PTC) in 
the framework of smooth solutions. Generalized 5-shock solutions resulting from Riemann problem 
are also investigated. Finally, using a kinetic scheme proposed in the literature in several areas, we 
validate such a numerical approach and propose a dedicated extension at the frontier of the moment 
space in the framework of both regular and singular solutions. This is a key issue for application 
fields where such an approach is extensively used. 

Key words. Quadrature-based moment methods; gas-particle flows; kinetic theory; particle 
trajectory crossing; entropic measure solution; frontier of the moment space 

subject classifications 76N15 (35L65 65M99 76M25 82C40). 

1. Introduction 

The physics of particles and droplets in a carrier gaseous flow field are described in 
many applications (fiuidized beds, spray dynamics, alumina particles in rocket boost- 
ers, . . . ) by a number density function (NDF) satisfying a kinetic equation introduced 



Crete numerical parcels of particles through a Lagrangian-Monte-Carlo approach or 
on a moment approach resulting in a Eulerian system of conservation laws on velocity 
moments eventually conditioned on size. In the latter case investigated in the present 
contribution, the main difficulty for particle flows with high Knudsen numbers (i.e. 
weakly collisional flows), where the velocity distribution can be very far from equilib- 



"Universite Paris Diderot-Paris 7 & Laboratoire J.-L. Lions, U.M.R. 7598 UMPC, Boite courricr 
187, 75252 Paris Cedex 05, France and Laboratoire EM2C - UPR CNRS 288, Ecole Centrale Paris, 
Grande Voie des Vignes, 92295 Chatenay-Malabry Cedex (chalons@math.jussieu.fr). This research 
was partially supported by an ANR Young Investigator Award (French Agcncc Nationale de la 
Recherche), contract ANR-08-JCJC-0132-01 - INTOCS - 2009-2013 

tLaboratoire EM2C - UPR CNRS 288, Ecole Centrale Paris, Grande Voie des Vignes, 92295 
Chatenay-Malabry Cedex and IFP Energies Nouvelles, (damienk@stanford.edu). D. Kah was funded 
during his Ph.D. Thesis by an IFPEN/EM2C CIFRE Ph.D. Grant. 

tLaboratoire EM2C - UPR CNRS 288, Ecole Centrale Paris, Grande Voie des Vignes, 92295 
Chatenay-Malabry Cedex (marc.massot@ecp.fr - Corresponding author. Present address : mmas- 
sot@stanford.edu - Center for Turbulence Research, Stanford University). This research was sup- 
ported by an ANR Young Inverstigator Award for M. Massot, ANR-05-JCJC-0013 - jeDYS - 2005- 
2009. We would like to thank the Center for Turbulence Research at Stanford University and its 
director, Parviz Moin, for their kind invitation during the Summer Program 2010 where the present 
work was completed and Laboratory EM2C for supporting the visit of C. Chalons at Stanford as 
well as Damien Kah's Master Thesis during which the study was initiated. 




Solving such a kinetic equation relies either on a sample of dis- 
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rium, is the closure of the convective transport at the macroscopic level. One way to 
proceed is to use quadrature-based moment methods where the higher-order moments 
required for closure are evaluated from the lower-order transported moments using 
quadratures in the form of a sum of Dirac delta functions in velocity phase space (see 
Yuan fc Fox (2010)] |Kah (2010)] and the references therein). 

Such an approach also allows for a well-behaved kinetic numerical 
scheme in the spirit of Bouchut |Bouchut fc al (2003)] ( See references from 



Laurent (2002) Freret fc al (2008) 



de Chaisemartin fc al (2008)] to 



de Chaisemartin (2009) 
Kah et al. (2010)] 



Massot fc al (2009) 



Freret et al. (2010) 



Yuan fc Fox (2010)] ) where the fluxes in a cell-centered finite-volume formula- 



tion are directly evaluated from the knowledge of the quadrature abscissas and 
weights with guaranteed realizability conditions and singularity treatment. Such 
a quadrature approach and the related numerical methods have been shown to be 
able to capture particle trajectory crossing (PTC) in a Direct Numerical Simulation 
(DNS) context, where the distribution in the exact kinetic equation remains at all 
times in the form of a sum of Dirac delta functions. Such methods can be extended 
to partially high-order numerical schemes ( |Vikas et al. (2010)] ). 

In another component of the literature devoted to multiphase semiclas- 
sical limits of the Schrodinger equation Jin fc Li (2003)] Gosse et al. (2003) 
Gosse fc Runborg (2005)] , the series of Wigner measures obtained from the Wigner 
transform for studying the semiclassical limits can be shown to converge towards 
a measure solution of the Liouville equation. Such an equation naturally un- 
folds the caustics and can generate the proper multiphase solutions globally in 
time. Two approaches have been used to solve this equation with a moment ap- 
proach, either the Heaviside closure |Brenier fc Corrias (1998)] as it is called in 
'Jin & Li (2003) Gosse et al. (2003")], o r, the one which is related to the present work, 



the delta closure (see Jin & Li (2003) Gosse et al. (2003")] and references therein). 



It leads to a weakly hyperbolic systems of conservation laws by taking moments of 
a Liouville equation exactly identical to the Williams-Boltzmann equation studies in 
gas-particle flows previously mentioned. Such approaches naturally degenerate to- 
wards the pressureless gas system of equation in the context of monokinetic velocity 
distributions |Massot fc al (2009)] |Kah et al. (2010)] |Kah (2010)] |Runborg (2000)] . 

Numerical algorithms in order to simulate such systems of conservations laws 
with the related delta closure or quadrature-based closure have been proposed in 
Jin & Li (2003) Gosse et al. (2003")] and [Desjardins et al. (2008)] independently, 



from the work for Bouchut & al (2003) Laurent (2002)] using naturally kinetic 
scheme with finite volume methods. However, many issues are still to be tackled 
in order to reach fully high order numerical schemes that preserve the vector of mo- 
ments inside or at the frontier of the moment space, thus leading to several possibil- 
ities of degeneration from a given number of abscissas to a lower one. In fact, such 
models are meant to capture a given level of complexity in the phase space which is 
fixed in advance by the number of moments and related quadrature nodes. In some 
particular situations, for perfectly controlled dynamics, it can be guaranteed that 
the solutions will remain smooth and consist in free boundary value (contact discon- 
tinuities) problem associated with switches between various numbers of quadrature 
abscissas. However, in most cases the numerical schemes have to tackle the possibility 
of singular solutions when the dynamics complexity goes beyond the one allowed by 
the model. In such cases the solution of the resulting system of PDEs is the viscosity 
solution and does not reproduce the exact dynamics in phase space and measure so- 



lutions are expected, for which we need a precise framework. More specifically, even 
if for the pressureless gas system, Bouchut (1994)] had set the correct mathematical 
background to define general entropic solutions, such a work had not yet been per- 
formed for higher order moment methods in the cited publications and no rigorous 
link has been provided between these and zero pressure gas dynamics. This is the 
purpose of the present contribution for both theoretical and numerical points of view. 

The paper is organized as follows. First we introduce the quadrature-based or delta 
closure velocity moment models for kinetic equations and focus on the four moment 
model in order to generalize to what can be done to higher orders but which would be 
difficult to expose due to algebra complications. The behavior at the frontier of the 
moment space is characterized as well as the mathematical structure of the system 
of conservation laws. We then define entropy conditions and provide, for smooth 
solution, the one-to-one kinetic-macroscopic relation. We then tackle the Riemann 
problem and define entropy measure solutions. Three examples of piecewise linear 
and singular solutions are then provided for which we rigorously identify the entropic 
character of the solution and which are then reproduced numerically. 

2. Quadrature-based velocity moment models for kinetic equations 

Consider the solution / = f(t,x,v) of the free transport kinetic equation 

d t f + vd x f = 0, t>0 7 xeR,veR 7 f(0,x,v) = f Q (x,v) (2.1) 

The exact solution is given by f(t,x,v) = f(0,x — vt,v) = fo(x — vt,v). Defining the 
i-order moment Mj = J v f(t,x,v)v l dv, i = l,...,N, N gN, the associated governing 
equations are easily obtained from (|2.1j) after multiplication by v % and integration 
over v, and write 

d t Mi + d x M i+1 =G, i>0. 

For the sake of simplicity, but without any restriction, we will focus our attention 
hereafter on the four-moment model 



( d t Mo- 
at Mi - 
d t M 2 - 
d t M 3 - 



d x M 2 = 0, 
d x M3 = 0, 
-d x M A = 0. 



(2.2) 



It will be convenient to write 



under the following abstract form 



BtM + B x F(M) = 0, 
with M=(Mo,Mi,M 2 ,M 3 )* and F(M) = (Mi,M 2 ,M 3 ,M4)*. 



(2.3) 



2.1. Quadrature inside the moment space 

This model is closed provided that M4 is defined as a function of M. In 
quadrature-based moment methods, the starting point to define this closure rela- 
tion consists in representing the velocity distribution of f(t,x,v) by a set of two Dirac 
delta functions, that is a two-node quadrature : 

f(t,x,v) = p\(t,x)8{v -vi{x,t)) + p 2 (t,x)S(v -v 2 (x,t)), (2.4) 

where the weights p\(t,x) >0, p 2 {t,x) >0 and the velocity abscissas vi(t,x), v 2 (t,x) 
are expected to be uniquely determined from the knowledge of M(x,i). Dropping 
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the (cc,i)-dependance to avoid cumbersome notations, such a function / has exact 
moments of order i = 0, ...,4 given by p\v\ + p 2 v\. The next step then naturally consists 
in setting 



M 4 = p\V x + p 2 v 2 



(2.5) 



where pi, p 2 and v\, v 2 are defined from M by the following nonlinear system : 

' M =p 1 +p 2 , 
M 1 =p 1 v 1 +p 2 v 2 , 
M 2 = piv{+ P2V2, 
M 3 =pivf +p 2 v%- 



(2.6) 



At last, it remains to prove that this system is well-posed, which is the mat- 



ter of the next proposition. We refer to Jin & Li (2003) Gosse et al. (2003) 
Desjardins et al. (2008)] for the proof. 



Proposition 2.1. System (2.3\) - (2.5\} - lt2.6\} is well-defined on the convex phase space 
n, also called the moment space, given by 

SI = {M = (M ,M 1 ,M 2 ,M 3 ) t ,M > 0,M M 2 - M\ > 0}. 

Moreover, setting U= (px,p2>piVi,p2V2) , the function U = U(M) is one-to-one and 
onto as soon as we set for instance vi>v 2 . Moreover we have < p\ < Mq and 
0<p 2 <M Q . 



Proposition 12.11 can be extended to the more general case of a 2fc-moment models, 
k>l. The velocity distribution is represented in this situation by a set of k Dirac 
delta functions, leading to M i = Y^ k j = iPjV 1 j, i — 0, ...,2k — 1, and M 2 k = J2 t j=iPj v j k ■ 

2.2. Hyperbolic structure inside the moment space 

The two-moment model, corresponding to k = 1 (one-node quadrature) writes 

d t p+d x pv = 0, 
d t pv + d x pv 2 =0, 

which is the well-known pressureless gas dynamics system. Recall that this model 
is weakly hyperbolic (the jacobian matrix is not diagonalizable) with v as unique 
eigenvalue, the characteristic field being linearly degenerate. Since there can be areas 
in the solution where a single quadrature node is sufficient at the frontier of the 
moment space in order to describe the dynamics, the solution in such zones will satisfy 
the previous system of two conservation laws. However we will first work inside the 
moment space and leave the behavior at the frontier for the next subsection. 

Actually, we will observe in the course of the next section that the four-moment 
model (|2.3[) is equivalent for smooth solutions (only) to two decoupled pressureless 
gas dynamics systems associated with (pi,piv\) and (p 2 ,p 2 v 2 ) respectively. Then 
(|2.3|) is expected to admit two eigenvalues v\ and v 2 and to be weakly hyperbolic 
with linearly degenerate characteristic fields, as stated in the following proposition. 



Proposition 2.2. ( f Jin & Li (2003%\Gosse et al. (2003)^ ) System $M$-$MM-WM 
is weakly hyperbolic on J! and admits the two eigenvalues v\ and v 2 , v\^v 2 . The 
associated characteristic fields are linearly degenerate. 
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Proof. For the sake of completeness, we propose here a direct proof of the eigen- 
values vi_ and v 2 of (|X5 ]) -(|21> ]) -P1) ]) , By (j2l)l) . we first easily get 



Mo=pl+P2, 
Ml—VlMo = P2(v2—Vl), 
Ml — V\M\ =P2V2(V2 — V\), 

M3-V1 M 2 = P2V%(v 2 -vx), 

and then, setting aa = viv 2 and a\ = — (vi +v 2 ), 

( Af o Mi \ (ao\ = _(M 2 
\M X M 2 ) y o~\ ) yM 3 

This system is invertible in the phase space SI (MoM 2 — Mf ^ 0) and uniquely defines 
ctq and (Ti with respect to M : 



1 



M M 2 



■M\ 



M1M3 
M1M2-M0M3 



(2.7) 



Then, we have 



= pivf + P2V2 

= pi vfvi +P2V2V2 

= {p\v\ +p 2 v 2 ){vi +v 2 ) ~ {pivf + p 2 vj)viv 2 



(2.8) 



= -M 2 a -M 3 a 1 , 

which finally gives M4 with respect to M. The Jacobian matrix J = VmF is given by 



J = 



/0 f 0\ 


I 


10 


with < 


f 


ya b c d J 





a = d Mo M i , 
b = d Ml Mj, 
c = d M2 M 4 , 
d = d M3 Mi. 



Using (|2.7[) and (|2 . 8[) . the calculations of the last row coefficients eventually lead to 



/ 1 

f 

1 

-2 ^-2 



\ 



\-o-g -2a Q a 1 -2a -o-f -2a x J 
Finally, the characteristic polynomial p(X) of J is easily shown to equal 

P (\) = (X-v 1 ) 2 (X-V2) 2 . 

This concludes the proof. □ 

Propositions |2~T1 and \2~2\ show that System (j2T3]) - (|23]) -([2~6 ]l is well-defined and 
weakly hyperbolic only on J7, which gives in particular V\ =/=V2 in the interior of the 
moment space. At a first sight, this might appear to be restrictive in the sense that 
one of the main objectives of the model is to allow particle trajectory crossing, that is 
in particular to deal with initial data consisting of two colliding particle packets such 
that V\ —V2 at each point initially (see for instance Section 7). Thus, in the last part 
of the present section, we characterize the behavior at the frontier T of the moment 
space when M o >0 : r = {M= (M , M 1 ,M 2 ,M 3 ) t ,M >0,M o M 2 -Mf = 0}. 
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2.3. Behavior at the frontier of the moment space 

As mentioned previously, it is rather natural to envision the coexistence, in a 
single smooth moment solution, of zones where the number of quadrature nodes are 
different. More specifically, we will examine the coexistence of zones where only one 
quadrature node is needed («i = v 2 ), that is where M = (Mo,Mi,M 2 ,M 3 )* with M > 
and M M 2 — M 2 = 0, and zones where M M 2 — M 2 > which are inside the moment 
space f2, whereas the vector of moments are smooth everywhere. 

After easy calculations in terms of p\, p 2 , v\ and v 2l the latter equality 
MqM 2 — Mf — writes p\p 2 {v\ — v 2 ) 2 —0, so that if the vector (pi,p 2 ,piVi,p 2 v 2 ) t 
exists, this actually corresponds to the case v:=vi — v 2 (still under the assumption 
pi 7^ and p 2 ^ 0) or to the case where one of the weights is zero. We also note that 
in both cases Mk = M$v k , whatever k in this case, so that the whole set of moments 
should be provided once Mp, and Mi are given, in close connection to the case of 
pressureless gas dynamics. There are in fact two possibilities : 

- either M = (Mo,M 1 ,M 2 ,M 3 ) t is such that M M 3 - M X M 2 ^ : in this case (|2l)j) 
cannot be solved and the vector (pi,p 2l p\Vi,p 2 v 2 ) t does not exist, 

- or M=(M ,M l ,M 2 ,M 3 ) t is also such that M M 3 -M 1 M 2 = : in this case (J2T6J) 
can be solved and we have v = v\ — v 2 — M\/Mq, together with p\ and p 2 defined 
by the one-parameter equation pi+ p 2 = Mq. As we will see just below, the choice 
pi = p 2 = Mo/2 is the most natural one when we have to deal with an isolated point 
at the frontier of the moment space. 

In order to justify the choice p\ = p 2 = M /2, we first observe that both condi- 
tions 

f M M 2 -M 1 2 = 0, 
lMoM 3 -MiM2=0, 

are equivalent to conditions 

(e = 0, (e = M M 2 -Ml 

\q = 0,' Wnere \ q = (M 3 M 2 - M?) - 3Mi (M M 2 - M 2 ) , 

and we consider p%, p 2 , v\ and v 2 as functions of (Mo,Mi,q,e) with Mq > 0, and e > C0. 
We then propose to study the asymptotic behavior of these functions when e— ^0 + , 
considering that M >0, M\ and q are fixed. Note that r = {(M ,Mi,e,g)*,M >0,e = 
0}. We get the following result. 



Lemma 2.3. Let be given M)>0 ; Mi and q. Then we have 

(M o ifq>0, ( 0ifq>0, 

lim p 2 = < 0ifq<0, lim pi={ M ifq<0, 



1 The definitions of e and q naturally comes out after noticing that setting p 1 = -jj^ , p 2 = jj^ ? 
vi = vi — -j^j , V2 = V2 — j^jj , solving 112.61 1 is equivalent to solving 

!l = Pl + P2^_ _ 
= Pi : ui +p 2 V2, 
q^-Lvf + Pzvl, 

with e = (MoM 2 -M 1 2 )/Af2 and q = ((M 3 M$ - Aff) - 3Mi (Mo M 2 - M\ ))/M$. 
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%£ifq>0, (+CX>ifq>0, 
lim v 2 = { -oo i/?<0, lim Wi = < # «/<?<0, 



MtifqX), ( Oifq>0, 

lim P2f2 =, \ Oi/g<0, lim /Oi^i=< Mi ifq<0, 

] ^ifq = 0, e ^° + l^tfg = 0. 

Proof. The admissible change of variables (Afo,Afi,Af2, M3) — > (Mo,Mi,e,q) al- 
lows to write after easy calculations 



Mi g+vV+ieS M e(i;iM -Afi) 

u i = i _ r H 777 ' /°2~ 



M 2M e ' P " V9 2 + 4e 3 ' 



and if q 7^ 



Mi o (l + si ff n(g)^l + 4e 3 /9 2 ) M e / UiM -Mi 
l 'i = 77" + T7— o > *°2 = I , 



Mq M e 2 g ysign(q) v /l + 4e 3 /g 2 

where we have set 



sign(q) = 



1 if q > 0, 
-1 if q<0. 



It is then an easy matter to get the expected results distinguishing between the three 
cases q < 0, q > and q = 0. It is then clear by a continuity argument that the proposed 
choice pi = p2 = Mq/2 when e = q = is actually natural. □ 

An important consequence of this lemma is that in the half plane e > 0, the region 
close to the frontier T for a non-zero q corresponds to abscissas going to infinity with 
arbitrary small weights. Moreover, when the velocity distributions at the kinetic level 
have compact support in the initial distribution, such a property will be preserved in 
the dynamics of the system and we want to be able to switch continuously from two- 
node to one-node quadrature without pathological behavior on abscissas and weights. 

Let us provide a first example where such a behavior is present. We consider a 
path in the moment space parametrized by the variable x, such that p\ = x 3 , P2 = 1, 
vi — l/x and V2 — 0. As x approaches zero, the smooth moment vector M=(l + 
x 3 ,x 2 ,x, 1)' has a very regular limit at the frontier of the moment space along the 
lines presented before with an unbounded abscissa. Indeed we have here e = x and 
q = 1 — x 3 approaches the fixed non-zero value of 1 . Note that if we replace the first 
weight by pi — x 4 , we still have an unbounded abscissa even if we converge toward 
the point (0,0) in the (e,q) plane (e = x 2 ,q = x(l — x 4 )). 

We thus have to find a framework in a subset of the plane (e,q) such that 
the limits are better behaved. A natural choice presented above is the line q = 
but it is too restrictive. In order to naturally introduce the relevant subset of 
f2, let us consider the other example with pi — ax 13 , p2 = l, vi—jx s and V2 — O, 
with a>0, 7>0, f3>0 and <5>0. As a; approaches zero, the moment vector 
M= (l + ax 13 ,ajx^ +s ,aj 2 x@ +2S ,Q7 3 i' 5+3 ' 5 )' reaches the frontier T of the moment 
spaco Two cases are interesting; firstly when (3 — 0, we reach the point (0,0) in 
the (e,q) plane asymptotically along the line q/e 3 ^ 2 = {l — a)/a 1 / 2 and no weight is 



2 The corresponding values of e and q are e = ofy 2 xP +2S and q = cry 3 ^ - ' -3 ' 5 (1 — ax^). 
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approaching zero, whereas the two abscissas are joining (see formulas in the proof 
above). Secondly, when <5 = 0, one of the weights is reaching zero, whereas the two 
abscissas remain different at a distance of 7 at the limit and we reach the point (0,0) 
in the (e,q) plane asymptotically along the line g/(Moe)=7 at the limit x — >0 (see 
again formulas in the proof above). We will prove in the following proposition that 
the proper framework is a symmetric cone in the (e,q) plane centered at the point 
(0,0) corresponding the |g/(M e)| <r/, where i] is a measure of the maximal distance 
allowed between the two abscissas. 



Definition 2.4 (Regular path). We define a regular path parametrized by x in 
the moment space, M x = (M 0x ,Mi x ,M2 X M 3x ) which admits a limit as x goes to zero 
and is at least C 1 up to the limit M . Moreover, we define its reduced second and third 
order moments e x — M 0x M 2x - M 2 X and q x = (M 3x M 0x -M? x )-3M lx (M 0x M 2x - 
M 2 X ). Its limit further satisfies e a — e x=o — and we assume e x>0 > 0, M Ox >v>0, 
\q x /(M Qx e x )\ <rj, where n>0. 

Proposition 2.5. We then have the following properties : 

• \im x ^ Q q x = <7o = 0. 

• the weights and abscissas admit limits p i0 = \mi x ^ p ix , v i0 = \im x ^ n v ix . If 
we assume that pio > for both i, then i>io = V20, or, if one weight approaches 
zero, such as pia = then we have \viq — v 2 q\ <f? and 77 is then a bound on 
the distance between the two abscissas. 

• M = (M 00 ,M 10 , M^o/Moo, Affo/M^o)*. 

Proof. It is first clear that lim x _>o q x = q = since \q x \ < nM 0x e x and e = 0. 
Then, easy calculations give 



so that denoting I = lim^o \q x \/(M 0x e x ) > 0, we clearly have v\ x — v 2x — > I when x — > 
and 77 represents an upper bound for v±o — V2o- Let us now distinguish between the 
cases I > and 1 = 0. Wc first note the following expression for qj (Moe) : 

-(uj 1 -uj 2 )(v 2 -v 1 ), u)i = Pi/M , uji+uj 2 = 1. 



M e 

If / > 0, one can write 

|P1x-P2*| = | I X 1 :XM 0x , 

M 0x e x \v lx -v 2x \ 

and this quantity clearly tends to Moo as x goes to zero. Which means that one weight 
approaches zero, /Oio = or p 2 o = 0, and the other one Moo (with p\ x +p 2x = Mq x ). 
If / = 0, it is clear by the following formula for vi (see the proof above) 



Mi <7+vV+4e 3 

ui = 1 

1 M 2M e 

that U10 = V20 = Mio/Moo- Using now the definition of p\ and p 2 (see again the proof 
above), one easily get 



p 2 _ \J 1 + 4e 3 / q 2 + sign(q) _ pi _ ^ 1 + 4e 3 / q 2 - sign(q) 

W2 ~Mo _ 2y/l+4e 3 /q 2 ' Wl ~Mo _ 2 v / l+4e 3 /q 2 : 



so that both weights have limits depending on the limit of q/e 3 ^ 2 . 

Clearly, M. = (M 0Q ,M 10 ,M? /M 00 ,M$ /M$ ) t , which completes the proof. □ 

A very important consequence of the previous proposition is the fact that along 
smooth paths inside the proposed cone which reach the point (0,0) in the (e,q) plane, 
the flux introduced in equation 12.31 is regular up to the frontier of the moment space, 
even if the mapping of M onto U is nolo 

Proposition 2.6. For any regular path in the moment space satisfying the assump- 
tions of the previous proposition, that is living in the proper cone in the (e,q) plane 
and reaching smoothly the point (0,0), the flux F(M) is continuous up to the frontier 
of the moment space and C 1 at (0,0) in any direction inside the proposed cone. 

Proof. The proof is rather straightforward when one has noticed the two equa- 
tions, the first of which is the expression (|2.8[) of Mj = — cr§M<x — 01 M3 as a function 
of M2, M3, <7o and <7i, and the second is the expression of ag and <Ti : 

/ g Mi (MA 2 e \ 

M eMo \M J M 2 

9 Mi 

- 2 — 

V M e M J 

Finally, following the same lines for the evaluation of the Jacobian matrix of the 
flux (see matrix J in section 2.2), it becomes clear that the flux is continuous and 
continuously diffcrentiablc in any direction inside the proposed cone. Besides, it can 
be easily seen that the expression of the flux as a function of (Mo,Mi,e,q)* becomes 

Wj_ ( q Mi (MA 2 e \ ( _e_ ( Mi\ 2 \ 
Mo \M eM \Mq) M$]\M$ \M J J 

\M e M J \M 3 ^ \Mj M M 2 j ' 

M4/M0 tends to (M10/M00) 4 when e goes to + which concludes the proof. □ 
Remark 2.1. Let us emphasize that in the various configurations we have proposed 
when the convergence toward the frontier of the moment space does not lie inside the 
cone \q x /{Mo x e x )\ <r\ in the (e,q) plane, the flux can dramatically loose regularity. It 
can either have a limit without being differ entiable or even not have a limit at all. The 
impact of the previous proposition thus becomes clear and sets the proper framework 
for solutions which will reach the frontier of the moment space. 

Remark 2.2. In the case Mq>0, at the frontier of the moment space within the 
previous proposed framework, we have e = 0, q — 0; the model is made of the two 
unknowns Mo and Mi and then degenerates to the usual pressureless gas dynamics 
which is weakly hyperbolic with a single eigenvalue v = M\/Mq. It should be noticed 
that for smooth solutions, the last two equations of system UTB on M2 — Ml/M^ and 
Ms = Mi/Mq are still satisfied with M4, = M±/Mq coherent with the previous limit 
obtained for the flux. As a consequence, we can notice, that at least for smooth 
solutions, the system of partial differential equations fH.ty) can describe the dynamics 
inside and at the frontier of the moment space. 

3 The mapping of M onto U will never be smooth at point (0,0) in the (e,q) plane since the 
limit of U will depend on the limit of q/(Moe), whereas the limit value of M in the proposed cone 
is always fixed. 
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3. Entropy conditions 

In this section, we will work in the interior of the moment space and we exhibit 
natural entropy inequalities for the following small viscosity system associated with 
(EH) : 

' d t M Q + d x M x =ed xx M , 
cWi + d x M 2 = ed xx M x , 
' d t M 2 + d x M 3 = ed xx M 2 , 
d t M 3 + d x Mi = sd xx M 3 , 

which gives in condensed form 

d t M + d x F(M)=ed xx M. (3.2) 
Throughout this section, we will consider smooth solutions only. We thus have 
d t M + 3d x M = ed xx M with J = V M F. 



Setting A = VuM, we then get 

dfU + A^JA^U^eA^a^A^U). (3.3) 
Our objective is to prove that 

dtV + d x q<ed xx r], (3.4) 

for a natural choice of couple (rj,q) given by 

f V = PiS(v 1 ) + p 2 S(v 2 ), , 35 x 
\ q = p x v 1 S(v x ) + p 2 v 2 S(v 2 ). 

Here S denotes a convex function from R to R, and we will especially consider the 
case where S(v)=v , a>0. Of course, the densities pi, p 2 and velocities v±, v 2 
involved in (13.51) are naturally defined by means of the one-to-one and onto function 
U = U(M). In the following and with a little abuse in the notations, we will consider 
without distinction r\ and q as functions of M or U. 



We first observe 



d t ri + d x q=d t r)(V) + d x q(U) 

= Vu^tU + Vug^U 

= V U 7/{-A- 1 JAa a: U + eA- 1 ^(A^U)} + VugaxU. 

The following two lemmas, the proofs of which are left to the reader, will be useful in 
order to estimate the entropy dissipation rate D defined by 

D = Vut^-A- 1 JAd x XJ + eA" 1 ^ {Ad x TJ)} + Vu^U. 

Lemma 3.1. The matrices J and A and A _1 JA are given by 

I 1 \ 

10 

1 

\-v\v\ 2v\v 2 (y\+v 2 ) -2viv 2 - (vi +v 2 ) 2 2(vi+v 2 ) J 
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A = 



1 10 0^ 
11 

—v\ —v\ 2v\ 2v 2 
\-2vf -2u| 3vf 3u| 7 



A JA = 



/ 


V o 



1 \ 
1 

2«i 
-v\ 2u 2 / 



Lemma 3.2. The gradients Vui] a?irf Vu9 are given by 



S(v 2 )-v 2 S (v 2 ) 
S'(vi) 
\ S'(v 2 ) J 



Vug: 



-viS\v 2 ) 
S{vi)+v x S («i) 

Vs(»i)+«iS'(»i)/ 



have 



V u g = V U 77A _1 JA. 



Before going on, let us make the following two remarks. We first note that thanks to 
the first lemma, (|3.3p with e = inside the moment space gives : 

' d t pi+d x p 1 vi=Q 1 
dtP2 + d x p 2 v 2 = 0, 
dtpivi - vfd x pi + 2v 1 d x p 1 vi = 0, 
t dtP2V 2 -v%d x p 2 + 2v 2 d x p 2 v 2 = 0, 



which is equivalent to 



d t pi + d x pxvi = 0, 
d t piv 1 +d x piv\=0, 
dtP2 + d x p 2 v 2 = 0, 
dtP2V 2 + d x p 2 vl=0. 



(3.6) 



Besides, at the frontier of the moment space, we obtain the pressureless gas dynamics 
on a single quadrature node. We then observe that for smooth solutions, thanks to 
the remark at the end of the previous section, the system (|2.3[) is nothing but either 
two decoupled or one single pressureless gas dynamics system of equations. We then 
observe that still with e = 0, D = 0, by lemma I3~2l in both cases : 

d t r) + d x q = 0, (3.7) 

that is for smooth solution both inside and at the frontier of the moment space. 
Let us go back to the case e > 0. We thus have the following equality, 

D = eV lJ r 1 A- 1 d x (Ad x U), 

from which it is natural to isolate ed xx n : 

D = eVu^A^cUAc^U) 

= eVu^a^U + e\7 v riA- 1 d x Ad x lJ 

= ed x {V\jnd x \5)-ed x {Vi]n)d x \J + e\7 v r] A' 1 d x Ad x \J 

= ed xx n + e(V v nA- 1 d x Ad x lJ -d x {\7ijri)d x lJ) . 



By lemma I3T2I giving Vu?7, we easily get 



d x (V v r ] )d x V = p 1 (d x v 1 ) 2 S"(v 1 )+p 2 (d x v 2 ) 2 S"(v 2 ). 
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It is now a matter to calculate Vu?7A 1 d x Ad x XJ — A t ('Vur]) t d x Ad x XJ . We first 
observe that 



cLAcLU = 



2pi(d x v 1 ) 2 + 2p 2 (d x v 2 ) 2 
\ 6pi vi (d x v 1 ) 2 + 6/9 2 w 2 (9^ v 2 ) 2 / 



so that only the last two components of A '(Vrjf))' are actually needed. Finally, easy 
calculations not reported here lead to 



A- t (V u r 1 ) t d x Ad x U = — 



1 



■{2p 1 {d x v 1 ) 2 X 1 + 2p 2 {d x v 2 ) 2 X 2 ), 



{ Vl -v 2 y 

where we have set 

X x = (vi - v 2 f (3(S(ui) - S(v 2 )) - 2{ Vl -v 2 )S' (yi) - (vt - v 2 )S' (v 2 )) , 
X 2 =-(«i- v 2 ) 2 (3(S(vi) - S(v 2 )) - ( Vl - v 2 )S' (vi) - 2(«i - v 2 )S' (%)) . 

The entropy inequality (I3.4[) is then valid if and only if 

VuyyA- 1 ^ Ad x U - d x (V v r))d x TJ < 0, 
that is, setting Si = S(vi), S^S (vi) and S i =S (v^), i = 1,2, 

(«i - « 2 ) 2 Pi (5^i) 2 (6(5! - 5 2 ) - 4(« x - v 2 )S' 1 ~ 2( Vl - v 2 )S' 2 + («! - « 2 ) 2 ^') 

+ 

(«i - « 2?P2{d x V2) 2 ( - 6(5! - S 2 ) + 2(«i - «2)^ + 4(«i - v 2 )S' 2 + ( Vl - v 2 ) 2 S 2 ) 



>0 



A sulhcient condition is given by 



6(5i - Sa) - 4(«i - u z)^; -2(«i - v 2 )5 2 + («i - f 2 ) 2 5i > 0, 
-6(5i - 5 2 ) +2(«i -«2)^ + 4(«i -u 2 )5 2 + (ui - v 2 ) 2 S 2 > 0. 



(3.8) 
(3.9) 



Let us focus for instance on the first inequality (the second one is treated in a similar 
way), and let us consider the left-hand side as a function of v 2 , for any given v\ : 

Fi{vi) =6(51 -S 2 ) -4(ui -v 2 )S 1 -2(«i -v 2 )S 2 + (vi - v 2 ) 2 5i . 
Differentiation yields 

J-;(« 2 ) = 4(5;-5;)-2( Wl -t; 2 )(5;'+5 2 '), 
^(«a) = 2(5i - 5 2 ) - 2(wi - v 2 )S 2 " , 
-7 r i'( u 2) = 2(w 2 -wi)5 2 '". 

It is then clear that Fi{vi)=J 7 1 (v\)=J 7 l (vi)=J 7 l («x) = 0. Then, provided that 
5 (w) > 0, Vi;, we easily get by a chain argument based on the sign of the derivative 
and the monotonicity property that .Fi(i; 2 ) >0, Vi>i,i; 2 . We have thus proved the fol- 
lowing proposition : 
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Proposition 3.3. Smooth solutions of VS. 2^ satisfy the entropy inequality \3.4]) 
for any entropy entropy-flux pair (r],q) defined by \3. 5\) provided that v—¥S{v) is a 
smooth function from R to R with nonnegative fourth- order derivative. In particular, 
the natural choice S(v) = v 2a with a > 2 is suitable. 

Remark 3.1. Any third- order polynomial may of course be added to the leading term 
of S , without changing the sign of the fourth-order derivative. However, if we focus 
on (strictly) convex functions v^S(v) in order to get a (strictly) convex entropy 
rj = ri(XJ), only first- order polynomials may be added without changing the convexity 
property. 

Remark 3.2. If we consider S(v) = l,v,v 2 ,v 3 , it is easily checked that VS. 9\) holds 
true with two equalities. In agreement with \S.1\) . these choices that lead to the pairs 
(r],q) — (Mi,Mi + i), i = 0,...,3, are admissible. 

Remark 3.3. In the case Mq >0, e = 0, q = with the additional conditions associated 
with the connection between the interior and the frontier of the moment space pre- 
sented in subsection \2.Sl we clearly have M^ — M^/Mq^ 1 . In this case, the entropy 
pairs clearly work also in such a case, for smooth solutions, and admits a smooth 
behavior at the frontier of the moment space in the cone we have defined previously. 

4. Kinetic-macroscopic relation for smooth solutions 

For smooth solutions, we established in the previous section that the four-moment 
model (|2.2j) is equivalent to the following two decoupled pressureless gas dynamics 
model 

d t pi + d x pivi = 0, 

d t pivi+d x pivf = 0, , 4 ^ 

d t p 2 + d x p 2 v 2 = 0, 
d t p 2 v 2 + d x p 2 v^ = 0, 



where p\, p 2 , v\ and v 2 are defined by the nonlinear system (|2.6j) . The aim of this 
section is to prove a rigorous equivalence result, still for smooth solutions, between 
this macroscopic model and the free transport kinetic formulation (|2.1I) when the 
velocity distribution is given by a set of two Dirac delta functions. This result 
is nothing but a generalization of the one given in Bouchut (1994)] for the usual 
pressureless gas dynamics model. 

Proposition 4.1. LetT>0 and pi(t,x), Vi(t,x) in C 1 (]0,T[xR) for i = l } 2. Let us 
define 

f(t,x,v)=^2pi(t,x)5(v-Vi(t,x)). 
i=l 

Then, pi and Vi solve Vt.2\) , or equivalently J^. lty , in ]0,T[xR if and only if 

d t f + vd x f = 0, in ]0,T[xRxM (4.2) 
in the distributional sense, i.e. if and only ifVc/) 6 C^°(]0,T[xR) and x £ C^°(M) 
i-t r 2 

^ Pl {t,x){d t <i>{t,x)+v t {t,x)d x <i){t,x))x{^{x,t))=Q. (4.3) 
i=i 



1.3 



Proof. Let us first assume that (I4.3[) holds true. Since the velocity functions 
are in particular locally bounded, one can successively choose \ £ C^°(M) such that 
x(v) = v k , k = 0, ...,3 for all v = Vi(t,x) and then get 




pj (t,x) (d t <f){t,x) + Vj(t, x)d x (p(t,x)) (v l (t,x)) k = 



for all 4> G C£°(]0,T[xR). Invoking the closure relation (|2.6[) . this gives the four- 
moment model (|2.2I) . or equivalently (|4.1I) . as pi(t,x) and Vi(t,x) are smooth functions. 
Conversely, let us assume that the partial differential equations of (|4.1[) are satisfied. 
Using the mass conservation equations, it is then usual to show that for i = 1,2 



Summing over £ = 1,2 and integrating past a test function cf> £ C£°(]0,T[xR) gives the 
expected result (|4.3p . This concludes the proof. □ 

This Proposition allows, as a corollary, to introduce a particular type of solution 
which will be denoted piecewise free boundary C 1 solutions. Such solutions correspond 
to a discontinuous connection from the interior of the moment space to the frontier 
through a contact discontinuity for which the Rankine Hugoniot solutions are trivially 
satisfied as well as the entropy conservation equation Q3.7p . 

Corollary 4.2. We consider the following distribution at the kinetic level f(t,x,v) = 
^2 i=1 Pi(t,x)6(v — Vi(t,x)) , where pi(t,x)>0 and V\{t,x) are taken as constants (or 
sufficiently smooth in some time interval), whereas p2(t,x) is zero except in a compact 
connected subset Kq at time t = o/R, where p2(0,x) > and V2(0,x) are two constants 
(or sufficiently smooth in some time interval) such that V2^v\. The resulting solution 
at the moment level exhibits two discontinuities at the frontier of the compact set K t 
which is the translation of set Kq at velocity V2; however the system \4-l\ as well as 
the entropy conservation equation \3. 9f ) are satisfied in the weak sense, that is the 
equations are satisfied in the usual sense where the solution is smooth and Rankine- 
Hugoniot conditions are satified at discontinuity points. 

Proof. Clearly, the moment solution will satisfy the system of conservation equa- 
tions everywhere except at the frontier of the K t set. The Rankine-Hugoniot jump 
conditions are trivially satisfied at the discontinuities where the mass flux associated 
to the first abscissa is pi(vi— V2) in the referential of the discontinuity and leads 
to zero jump conditions for the part of the flux associated to the first abscissa by 
continuity, whereas the mass flux associated to the second abscissa is zero, as usual 
in contact discontinuities, which also allows to conclude. The same path allows to 
conclude that for any entropy-flux pair, the conservation equation is satisfied in the 
weak sense. □ 

Let us underline the fact that in the region where the second weight is zero outside 
the compact set K t , we have used so far the convention that in such a region where a 
single quadrature node is to be found, the two weights are equal and the two abscissas 
are equal. The results proposed in the previous corollary are of course independent of 
such a choice since the point at the frontier of the moment space is isolated. Besides, 
such a corollary can be extended to as many quadrature nodes as needed as long as 



Pi(dtVi+Vid x Vi) = 0, 



and then multiplying by x' 



' for any smooth function x, 



dtPix( v i)+d x piX{ v i) v i =0. 
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the number of nodes allows to describe the dynamics at the kinetic level at any point 
and time. Finally, the collision of two particle packets presented in subsection 16.11 
satisfies the assumptions of Corollary 14.21 and will be an entropic solution. 

5. Riemann problems and entropic measure solutions 

In this section, we focus on the Riemann problem, which is associated with the 
inital condition for two constant states Ml and in Q. 

M(x ' 0) = \M^ifx>0, (5 - 1} 



The solution of 



5.1[) is sought in the form 



'Mi tfx<a L t, 
M{(t)d{x-a L t) if x = a L t, 
M(x,i)=^M+ iia L t<x<a R t, 
M s R (t)5(x-a R t) if x = a R t, 
M R ifx> a R t, 



(5.2) 



which corresponds to the juxtaposition of two Dirac delta functions with mass Mj| 
and position x — apt, f3 = L,R, and separated by a constant state M* in Q. Here, 
ap denotes a real number, we have mp(t) >0, m,3(0) = and (3 = L 7 R, and M^(t) is 
defined by 



M 5 Jt): 



fmp(t) \ 
m (j {t)op 
mp(t)a 2 f 



(5.3) 



We introduce the following natural definitions of (entropy) measure solutions. 

Definition 5.1. Let n^ piS(v 1 ) + p 2 S(v 2 ) and q = piviS(vi) + p 2 v 2 S(v 2 ) with 
S(v)=v 2a , aSN. We say that \5.2]) is a measure solution of \2.3]) - ^5~l]) if and 
only if d t r] + d x q = in V (]0,oc[xR) for S(v) — v 2a with a = 0,i,l,§, that is for 
(r),q) = (Mi,M i+ i), i = 0,l,2,3. We say that \5.2\) is an entropic measure solution 
of pj [) -fib ]) if and only if d t r) + d x q<0 in X>'(]0,oo[xM) for S(v) = v 2a with 
a = Q,i,l,5 anda>2. 

We can prove the following equivalence result. 

Theorem 5.2. The solution given by 115. 2\) is a measure solution of i2.3\) - ^S~l]) iff 
a L (M L - M*)t - (F(Mi) - F(M*))i + M| (t) = 

a R (M, - M R )t - (F(M*) -F(M R ))t + M s B (t) = 



(5.4) 



The solution given by (5.2]) is a entropic measure solution of i2.3\) - ![5~l)) if and only 
if in addition, for all S(v) = v 2a with a > 2 : 



c7 L (r ? (M i )-r 7 (M Jt ))t-(g(M i )-g(M*))t+r ? (Mi(t))<0 

a R ( V (M,)- V (M R ))t- (q(M*)-q(M R ))t + r)(M R (t)) < 0, 
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(5.5) 



Proof. We first recall that dtr] + d x q = in T> (]0,oo[xR) means that for any 
smooth function with compact support <p £ C£°(]0,oo[xR), we have 

<d t ri + d x q,ip>=-<ri,dtip> -<q,d x tp>= 

~ Io° J-™ v( x ,t)dtf(x,t)dxdt ~ J °° f+™ q(x,t)d x ip(x,t)dxdt = 0, 



where by definition (|5.2 



rj(x,t) = < 



x < a L t, 

r]{M{(t))S(x-a L t), x = a L t, 



( V (M L 
{ V(M R ) 



(q(M L ), 



x < a^t, 



V (M R {t))d(x-a R t), x = a R t, 
x > a R t, 



°L<J< o R , q(x,t) = < g(M*), 



q(M{(t))5(x-a L t), x = a L t, 

0~L<T <VR, 



q(M 5 R (t))S(x-a R t), x = a R t, 



U(Mb), 



x > a R t. 



Here we note that q(M s p (t)) = o-pr)(M S p(t)) with j 1 (M 6 fj (t))=m [j {t)S{ap), /3 = L,R. For 
all <peC c °°(]0,oc[xR) we thus have 



<d t r] + d x q,(p> = 



dt 



J -oo 

oo rcrRt 



dt 



Ja L t 
oo poo 



r)(M L )d t cp(x,t)dx- / rj{M d L (t))d t ip(a L t,t)dt 
Jo 

poo 

r){M*)d t iptx,t)dx- / r](M R (t))dM<T R t,t)dt 
Jo 



oo ra L t 

dt I ri(M R )dt<p(x,t)dx- I dt q(M L )d x p(x,t)dx 

n JlTRt JO J-oo 

oo poo r&Rt 

q{M{(t))d x ip(a L t,t)dt- / dt q(M*)d x (p(x,t)dx 



o 



art 

oo poo poo 

6 





oo 

q{M d R (t))d x ip(a R t,t)dt- I dt I q(M R )d x <p(x,t)dx 

n JO JdRt 



that is, using in particular q(M.0(t)) = <r^?7(M^(i)), 



<d t r] + d x q 7 (p>^ 



f°° d f°° 

dt V {M L )(— tf(x,t)dx-a L <p(t,a L t))- q{M L )<p{cr L t,t)dt 

JO O^J-oo Jo 

poo ^ r&Rt 

J dtri(M+)( — J (p(x,t)dx - a R ip(t,cr R t) + o- L ip(t,a L t)) 



q(M*)(v(o R t,t)-y(<j L t,t))dt 



dt V (M R )( 



(p(x,t)dx + a R tp(t,a R t)) 



q{M R )<p(<r R t,t)dt 



d 1 r°° d 

V {M s L {t))-[ V (<j L t,t)]dt- I V (M s R (t))-[ V (a R t,t)]dt 



dV 



16 



rOO /'OC 

= / <p(a L t,t)(<T L r)(M L )-q(M L ))dt+ / ^(<T R t,t)(a R r 1 (M+)-q(M+))dt 

J 9oo i%o 

- \ V {<7 L t,t){a L n{M.*)-q{M*))dt- \ ^(a R t,t)(a R n{M R )~q(M R ))dt 
Jo Jo 

poo j poo , 

-J V<Wi(t))#[<P(<rLt,t)]dt- J V (M 5 R (t))-y(a R t,t)]dt 

(p(a L t,t)(a L (i](M L )-ri(M ir )) - {q(M L ) - g(M*)))# 
p(0*t,t) (o*fa(M*) - »7(Mk)) - (?(M*) - g(M K )))dt 

) T /"OO 7 

v (M s L (t))-M<j L t,t)]dt- / ^(M^(t))^[ V (<7 H t,t)]dt 



/•OO i 

= / ^( C r i t,i)-( ( r i (r 7 (M i )-r ? (M*))i-(5(M i )-g(M i ))t + ry(M£(t)))df 

poo 

+ jf ^N,t) s (<r B WM 1 )- f) (M fi ))f-( !/ (M 1 )- 9 (M JJ ))H!)(Ml(t)))& 

Then, since 771/3(0) =0 and then ^(M^(0)) = 0, f3 = L,R, it is clear that ([5^2]) is a 
measure solution of (|2.3p - (|5.1|) if and only if (|5.4p is valid for all t > 0, and an entropic 
measure solution if and only if in addition (|5.5[) holds true for all t > and all 7/ = 
piS(vi) + p 2 S(v 2 ) and g = /Oi^iS^i) +p2V2S(v 2 ) with S(tj) = tj 2q , a>2. □ 
Remark 5.1. Let us recall that r](M. S g(i)) = mp(t)S(o~p). Since i f 5. 5\) is made of 
equalities when S(v) = 1 and S(v) =v (we get in these cases the first two components 
of \5.4\ l), the validity of i5.5\) for all S(v) —v 2a , a>2, is equivalent to the validity of 

o-L (i)(Mi) - rj(M*)) - (q(M L ) - g(M*)) < 

(5.6) 

a^M*) - V (M R )) - (g(M*) -g(M R )) < 0, 

/or aiZ = t; 2 " - ^ - a£ ^ , a > 2 . 
6. Examples of entropic solutions 

In this section, we propose three particular entropic solutions. The first one 
models the collision of two particles packets with free boundary C 1 smooth solution, 
that is a solution for which an exact link with the kinetic level is preserved and 
for which the entropy equation is exactly satisfied. In such a situation the four- 
moment model does not develop 5-shock Dirac delta functions and is actually able to 
properly represent the crossing of the two packets which correspond to the dynamics 
at the kinetic level. The second one models the collision of four particles packets. In 
this case and as expected since the number of moments is set to four, the entropic 
solutions involves two <5-shock Dirac delta functions singularities. Whereas the first 
case corresponds to a connection from the interior of the moment space to the frontier 
through a contact discontinuity, or free boundary solution, resulting in an isolated 
point at the frontier, we consider in a third example a smooth connection to the 
frontier of the moment space, such that the point at the frontier is an accumulation 
point of a trajectory inside the moment space. 

6.1. Collision of two particles packets 

We consider a Riemann initial data (I5.ip where Ml = M(Ul) and M# = M(Ur) 
are such that 
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Ui =2 



/ PL \ 

PL 

PLVL 

\plv l J 



and 



V R = - 



( P* \ 

PR 

PrVr 
\prVrJ 



for two given densities pl>0 and pr>0 and velocities vl > and vr < 0. We recall 
that the function M = M(U) is defined by ([2~6]) . We define 



M L if x < v R t, 
M(x,t) = { M* if v R t<x<v L t, 
M R iix>v L t, 



3.1) 



with M* = M(U*) given by U* = (pltPRiPlVltPrvr) 1 ■ Our objective here is to prove 
that the following Dirac delta functions free solution is an entropy solution of (|2.3j) - 
(|5.ip . Conditions (J5T3J) and (|53j) write here 



5.2) 



and 



j v R (M L - M*) - (F(M L ) - F(M*)) = 
1 v L (M* - Mr) - (F(M*) - P(Mh)) = 0. 

Ui? (t7(M l ) - rj(M*)) - (?(Mx) - ?(M*)) < 
«i (»/(M*) - J7(M fl )) - ( g (M*) - q(M R )) < 0, 



(6.3) 



with rj(\J) — piS(vi) + p2S(v2) and g(U) = piUiS'(wi) + p2V2S(v2) for all S^u^ir" 
with a > 2. We will focus only on the first equality of (|6.2p and the first inequal- 
ity of (|6.3p . the second ones being treated similarly. We clearly have 
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It is then clear that the first equality of (16.21) holds true. Let us now check that the 
proposed Riemann solution fulfills the entropy condition. We clearly have 

v R (r,(M L ) -r?(M*)) - (q(M L ) -?(M*)) 
= vr(plS(v l ) - (plS(v l ) + prS(v r ))) - (p L v L S(v L ) - (plVlS(v l ) + Prv r S(v r ))) 

= 

which allows to prove that the proposed solution is an entropic smooth solution. 
6.2. Collision of four particles packets 

We consider a Riemann initial data fl5.1[) where Mi = M(Ul) and M^ = M(Ur) 
are such that 



( P \ 
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\PV2 j 



and 
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for a given density p > and two velocities v 2 > V\ > 0. We have 
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We define 



M(x,f)= < 



( M L itx<-at, 
M s L (t)5(x + at) i£x=-at, 
M* if -<rt<x<at, 

M 5 R (t)8(x-at) iix = at, 



(6.4) 



M; 



if x > at, 



with a > 0, M* = M(U*) given by 



U *=2 



and M s L (t), M 5 R (t) given by 
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with m(t) >0. The generalized Rankine-Hugoniot jump conditions (|5.4[) write here 
-a(M L -M t )t-(F(M i )-F(M,))HM[(t)=0 

rr(M*-M Jl )t-(P(M*)-P(Mfl))t+M^(t) = 0. 
that is, equivalcntly 



2aM* - (t(Mj,+Mj{) + (F(M fl ) - F(M L )) - 
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V o J 



2F{M,)+a(M R - M L ) - (F(M L ) + F(M L )) - 



m(t) 
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2cr 





(6.5) 
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This is made of eight equalities, four are trivial (zero equals zero), so that four are 
left to determine the four unknowns p+, u*, a and m(t)/t. We propose below to 
numerically solve this nonlinear system for a specific set of values for p, v\ and v 2 . 

Remark 6.1. We conjecture existence and uniqueness of a solution to this nonlinear 
system. Indeed, the initial condition involves four different velocities while the model 
is able to represent two different velocities only. More precisely and by analogy with the 
usual pressureless gas dynamics model (see \Bouchut & al (2003)^ , fBouchut (1994^ ), 
velocity v\ is to "bump" into velocities —v± and —V2 to create a first Dirac delta 
function. By symmetry, another Dirac delta function is expected. 

Regarding the entropy inequalities (|5.5j) . we first remark that for S(v) —v 2a , 

r,(M L ) = r 1 (M R ) = ^(S(v 1 ) + S(v 2 )) 
i7(M*) = /j*0(v*), 
r,(Mi(t))=T,(M s R (t))=m(t)S(a), 

while 

q(M L ) = -q(M R ) = ^(v 1 S(v 1 )+v 2 S(v 2 )), g(M*) =0. 

As an immediate consequence, both inequalities in (|5.5[) are equivalent and the entropy 
condition is 

- § (S(vi) + S(v 2 ))) -E( Vl s( Vl )+v 2 S(v 2 )) + ^S(a) < 0. (6.6) 

A concrete example. We propose to take p = l, «i =0.8 and v 2 — 1.2. Numerically 
solving (E31) gives = 1.88265, = 1.06026, a = 0.87983 and m(t) /t = 0.22342. If the 
left-hand side of (|6.6[) . which represents the entropy dissipation rate associated with 
S(v) = v 2a , is denoted D(a), a simple calculation gives for instance D(2) = —0.27324, 
D(3) = -0.86854, D(4) = -1.88678,... The proposed solution (JO]) is then actually an 
entropy measure solution of the four-moment model. Such an exact entropic solution 
will be used in the following to prove the relevance of the numerical scheme proposed 
hereafter with respect to the exact solution when singularities are present. 

6.3. Piecewise linear solution connected with the frontier of the mo- 
ment space 

As a last example, we introduce a piecewise linear solution which allows to connect 
zones inside the moment space with zone at the frontier within the proper framework 
introduced in subsection 12.31 Particles are initiated in the domain [0,0.5]. In the 
domain [0,0.1], a monomodal velocity distribution is reconstructed, with vi=v 2 = l, 
and p\ — p 2 = 0.5. On the contrary, a bimodal velocity distribution is reconstructed 
in the domain [0.1,0.4], with p\= p 2 = 0.5 and for the abscissas, Vi = 1+ I Q °' 1 and 
v 2 = 1. The initial conditions are represented in Fig. 16.11 

The ground difference with the first test case is that the transition between the 
two zones is smooth, and so the numerical strategy to account for this transition is 
important. The analytical solution of this problem, in smooth areas, consists of a 
decoupled transport of each of the quadrature nodes as two indcpcndant pressureless 
gas as showed in system (|3.6p . This comes from the fact that the number of Dirac 
delta functions reconstructed from the moments, two in this case, is always sufficient to 
capture the problem dynamics. The equivalence between the kinetic and macroscopic 
equations is preserved, and the solution in terms of moments satisfies the entropy 
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Density 




Fig . 6.1. Moment dynamics for a free boundary connecting areas where e = and e > . Initial 
conditions. Top-left: Mo. Top-right: M\, Bottom-left: weights. Bottom-right: abscissas. The solid 
line corresponds to (pi,vi), the dashed line with circles to (p2,V2). 



equation. Therefore, the solution is the superposition of a translation at constant 
velocity v% = 1, and a transport with the following velocity field: 

1 i£ [0,0.1] 
vi={ 1 + ^ir x e[0.1,0.4] (6.7) 
2' xe 0.4,0.5 



The analytical solution, displayed in Fig. 16.21 has two fronts at x = 0.2 and x = 0.9 
moving with a velocity v = 1 and v = 2 respectively. The first front corresponds to 
particles initiated with velocity v = l between x = and x = 0.1. The second front 
corresponds to particles initiated with velocity v = 2 between x = 0.4 and x = 0.5. The 
square wave between x = 0.8 and x = 0.9 is the final location of the particles initiated 
with velocity v = 2 between x = 0.4 and x = 0.5. The value of p\ between x = 0.3 and 
x = 0.8 corresponds to the expansion of the density field due to transport with a linear 
velocity field with positive slope. The value of the density in that area is the solution 
of the equation dtpi +v\d x p\ = —pid x vi, which yields pi =0.3 at time t = 0.20. 

According to the initial conditions, both weights have the same profile, the quan- 
tities q/(Moe) and qje 3 ! 2 are null. At time t = 0.2, the weight profiles are different in 
the interval [0.3,0.8] corresponding to the expansion of p\. Therefore, q/e 3 ^ 2 is non 
null, as well as jg-^, since the velocities have also different values, and can be exactly 
calculated (see Fig. 17.101) . 

7. Numerical simulations via kinetic schemes 

This section is devoted to the discretization of (|2.3[) - (12.5[) - (|2.6p . As already stated, 
we use as a building block a natural first-order kinetic scheme already proposed in the 
literature |Jin fc Li (2003)[ |Gosse et al. (2003)[ |Desjardins et al. (2008)] and briefly 



recalled here for the sake of completeness. 



4 In the interval [0.8,0.9], although p2 should be null (the square wave at velocity V2 = 1 should be 
bounded between the front a; = 0.2 and a; = 0.7) we computed p2 = 0.5 and V2 =2. This is consistent 
with the conditions in Section 12.31 for moment vectors at the frontier of the moment space. 
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Fig. 6.2. Moment dynamics for a free boundary connecting areas where e = and e > 0. Solution 
at t = 0.2. Top-left: Mo. Top-right: M\. Bottom-left: weights, Bottom-right:abscissas. The solid 
line corresponds to (pi,vi), the dashed line with circles to (p2,V2)- 



Let us first introduce a time step At>0 and a space step Ax>0 that we 
assume to be constant for simplicity. We set X = At/Ax and define the mesh 
interfaces Xj + \/ 2 =jAx for j e Z, and the intermediate times t n = nAt for n€N. 
In the sequel, M" denotes the approximate value of M at time t n and on the cell 
C j = [x j _ 1/2 ,x j+1/2 ). For n = 0, wc set M° = ^ / * ^ M Q (x)dx, jeZ, where 
Mo (a;) is the initial condition. 

Let us now assume as given (M™) je z the sequence in Q of approximate values at 
time t n . In order to advance it to the next time level t n+1 , the kinetic scheme is 
decomposed into two steps. 

First step : transport (t n — >t n+1 ~) 

We first set U" = U(M") and define the function (x,v) — >f n (x,v) by 

f n (x,v) = (pi)]S(v-( Vl )]) + (P2)]5(v-(v 2 ^), V (x,v) GCj xR, j G Z. 
We then solve the transport equation 

jd t f + vd x f = 0, (i^jelxM, 

\f(t = 0,x,v) = f n (x,v), [ '- L> 

the solution of which is given by f(t,x,v) = f n (x — vt,v). 
At last, we set f n+1 -(x,v) = f n (x-vAt,v). 

Second step : collapse (t n+1 ~ -> t n+1 ) 

The first four moments at time t n+1 are now naturally defined by setting 

i r x i+i/2 r+oo 
(M^ +1 = — / v\f n+1 -(x,v)dvdx. 

L± X Jx j _ 1/2 J -OO 
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Then, we have M] +1 = ((M )J +1 ,(M 1 )^ +1 ,(M 2 )] +1 ,(M 3 )] +1 ) t for all j eZ, which 
completes the algorithm description. 

Remark 7.1. It is easy to see that this scheme preserves the moment space ft, see 
for instance fDesjardins et al. (2008 fj . 

Remark 7.2. Under the natural CFL condition At max jeZ ((ui)™, (v 2 )") < CFL Ax, 
with CFL<\, integrating |7. lty over (t,x,v) G (0, At) x Cj x K and against v l , i = 
0,...,3 easily leads to the equivalent update formula 

Mr 1 =M--^(F^ 1/2 -F»_ 1/2 ), jeZ, 

where we have set F» +1/a = ((Mi)? +1/a> (M a )» +1/2 , (M 3 )? +1/2 , (M 4 )? +1/2 )' and 
(Mi)? +1/2 - (M<)?+ /a + (M 4 )^ 1/2; and wjft 

(M^ 1/2 = (pi)? +1 min(0>^^ 

( M 0;i /2 = (Pi)>ax(0,K)^ 1 )(K) J ") I - 1 + (p 2 )>ax(0,( W2 )^ 1 )(( W2 );0 1 " 1 . 



7.1. Numerical quadrature strategy at the frontier of the moment space 

This paragraph addresses the issue of how to numerically handle the transition 
between a vector in the interior of the moment space, and a vector lying at its frontier. 
For an isolated point at the frontier of the moment space T, there is no specific problem 
since we use a single quadrature node and the quadrature is not a problem. The two 
problems we have to face are related to the preservation of the cone in which we 
envision to work in section 2 as well as to deal with finite precision algebra in the 
neighborhood of the point (0,0) in the (e,q) plane. 

Consequently we introduce two constants in the numerical quadrature we use. 
First, for finite precision algebra and in order to avoid numerical errors, we only 
evaluate the two quadrature nodes when e/M^ > ei, where t\ is a small number related 
to machine precision. Under this threshold, the velocity dispersion is considered null 
(e = 0), and the set of Dirac delta functions are reconstructed as suggested in section 
2: pi=p 2 = A/o/2, vi=v 2 = Mi/M . 

Secondly, since we want to deal with compactly supported velocity distribution 
at the kinetic level, we will introduce another constant, rj which is a bound for the 
distance between the two abscissas. In fact we would like to impose to the solution to 
remain inside the cone in the (e,g) plane : < 77. It will be shown in the following 
examples that the cone is in fact automatically preserved by the proposed algorithm 
and that such a bound does not have to be imposed but is satisfied by the numerical 
solution. 

If -jJ|L > ? y j then we set q/(M e) — sign(q)rj, so that : 

~-^ = sign(q)T], q = sign(q)nM e, M 3 = q+M 1 M 2 /M . 
M e 

As illustrated in Fig 17.11 and explained in Section 12. 31 the variable q naturally lives 
in a cone delimited by the straight lines of slope ±77 and becomes null when e < t\. 
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Fig. 7.1. Handling of the moment space border with the admissible cone. 



Let us remark that limiting does not lead to a limitation on the quantity 
-§j2 in such a way as to allow the possibility to reach large ratio in terms of weights 
(one weight can approach a zero value while the quantity -^pz grows indefinitely) 
but without allowing the abscissas distance to grow beyond of fixed limit naturally 
inherited from the initial solutions at the kinetic level. We will come back to this 
point in the following result section. For the simulations we present, the parameters 
are such that: t\ — 10 -9 and r; = 2. 

7.2. Numerical results 

This section is devoted to numerical illustrations of the two Riemann problems 
and one dedicated to the case of a free boundary connecting zone where e = and a 
zone where e > 0, discussed in Section 6. 

In all the figures, we choose to represent p\ and v\ by solid lines, and p2 and vi 
by lines with circle markers. In the representation of weights and abscissas, values 
have to be assigned for v\ and vi : we thus decide that v\ is the maximum of the 
relative values of velocity. 

Two packet collision. Figure 17.21 represents the initial conditions for the 
two particle packet case. Figures 17.31 and 17.41 presents the numerical and analytical 
solutions for the two particle packet case with p^ = p^ = 1, and Vt, = 1, vr — — 1. The 
computation is run with a 1000 cell grid on the spatial domain [0,1], with CFL = 1. 
The length of each packet is 0.4 (pi = p2 = v\ = v-i = for x < 0.1 and x > 0.9) and the 
two packets start to collide exactly at time t = 0. Moving in opposite direction one 
across the other, with the same opposite speed, they then overlay and we note in 
particular that pi — p2 = 1 and v\ = V2 = in the mixing zone (see for instance the 
plots at time £ = 0.1). As expected, they finally get separated again and we note that 
a perfect agreement is obtained with the exact entropic solution. 

Four packet collision. Figures 17.51 presents the initial conditions. Figures 17.61 
and 17.81 present the numerical and analytical solutions respectively for the moments 
and the weights and abscissas. The computation is run with a 1000 cell grid on 
the spatial domain [0,1], at CFL = 1. Here, we observe the presence of two Dirac 
delta functions as already discussed in Section |6] The agreement between exact and 
numerical solutions for the moments is very good, showing that the numerical solution 
converges to the analytical one. The disparities encountered between the analytical 
and numerical solution in the case of the weights and abscissas are due to the fact that 
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Fig. 7.2. Initial fields of weights (left) and abscissas (right) for the two particle packet case. 
The first quadrature node is represented by solid lines whereas the second node is represented by 
circles. 
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Fig. 7.3. Results for the two packet case at time t = 0.1. Top: Numerical results for Mo (solid 
line) and Mi (dashed line) (left) and for M2 (solid line) and M3 (dashed line) (right). Middle: 
Analytical result for weights (left) and abscissas (right). Bottom: Numerical results for weights 
(left) and abscissas (right). The solid line corresponds to the higher abscissa, the dashed line with 
circle to the lower one. 
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Fig. 7.4. Results for the two packet case at time t = 0A. Top: Analytical result for weights 
(left) and abscissas (right). Bottom: Numerical results for weights (left) and abscissas (right). The 
solid line corresponds to the higher abscissa, the dashed line with circle to the lower one. 



the mapping U(M) is discontinuous at the moment space border. In area of numerical 
diffusion (xw0.22 and x«0.88), when the dispersion e is under the threshold ej., 
weights and abscissas are reconstructed as explained in section 2. 

The wave propagating at velocity 1.2 and separating the constant states (vi — 
0.8,i;2 = 0) and (vi — 1.2, — 0.8) is steep and coincide with the analytical wave 
whereas the wave propagating at velocity 0.8 is actually smooth since the CFL number 
is based on the highest value of velocity, which is 1.2 in this studied case. The same 
explanation holds for the symmetric jump at location x — 0.78. Meanwhile, because of 
the conservation of the velocity moments, the numerical velocity jump (at x = 0.154) 
happens before the analytical velocity jump (at x = 0.18). One can here notice that 
the quadrature method provides the expected value of velocity in the numerical dif- 
fusion zones. The same explanation holds for the different velocity jump locations 
between the analytical and numerical solution at a; = 0.82 and a; = 0.845. The same 
phenomenon is responsible for the disparities between the analytical and numerical 
solutions at the (5-shocks locations, at x = 0.4 and a; = 0.6. 

Figure 17.71 displays the final profile of the quantities and Thus, 
has significative values in areas where the abscissa distance as well as the weight 
difference are important, whereas reaches high value in areas where the weight 
ratio is important. Since in the domain, except for the singularities, the weights have 
the same value, both quantities are equal to zero. At the singularities, -Ay is roughly 
proportional to the square root of the weight ratio, and jg-g is bounded, accounting 
for the fact that the velocity field is bounded. 

We have thus provided numerical simulations in the two cases for which we 
have at our disposal an analytical entropic solution, either in the piecewise contant 
case, or in the singular case where <5-shock mesure solutions are present. In the 
first case, the crossing of the two droplets monokinetic packets is very properly 
reproduced without numerical diffusion since we work at CFL one, even if this 
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Fig. 7.5. Four packet case with PL = PR = 1, o,nd vi = l.2, 1)2=0.8. Initial conditions. Top: 
Mq (left) and M\ (right). Bottom: weights (left) and abscissas (right). The solid line corresponds 
to the set (pi,v\), and the dashed line with circles to the set (p2,V2). 




Fig. 7.6. Four packet case. Results at t = 0.1. Analytical (solid line) and numerical (dashed 
line with circles) solutions. Top: Mq (left) and Mi (right). Bottom: M2 (left) and M3 (right). 



is not symptomatic of the numerical diffusion such methods will encounter with 
a first order method in realistic configurations Kah el al. (2010)] . In the second 
case, the numerical method is able to capture the creation of the measure singular 
solutions associated to the fact the we have limited the number of quadrature node 
to two. With this node number limitation, the proper physical solution, in the 
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Fig. 7.7. Four particle packet case, at time t = 0.1. Left: ratio q/(Moe). Right: ratio q/e 3 / 2 . 
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Fig. 7.8. Four packet case. Results at 4 = 0.1. Top: Analytical weights (left) and abscissas 
(right). Bottom: Numerical weights (left) and abscissas (right). The solid line corresponds to the 
set (pi,vi), and the dashed line with circles to the set (p2,V2). 



infinite Knudsen number limit, where the various droplet packets cross without 
interacting, differs from the entropic solution of the system of partial differen- 
tial equations (|2.2[) obtained through the quadrature-based closure. This is the 
same type of behavior as seen in the case of pressureless gas dynamics at a lower level. 



Free boundary case. The last test case, explained in subsection 16. 3[ assesses 
the ability of the method to solve free boundary cases connecting in a continuous 
manner states lying in the interior of the moment space and at its frontier. The chosen 
grid contains 400 cells, CFL number is set to 0.98. This value of the CFL number is 
chosen is order to prevent high frequency instabilities to occur. The computation is 
run until t — 0.2. The analytical solution has been provided in subsection 16.31 

Results are displayed in Fig. 17.91 Let us hrst focus on the fronts present at 
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a; = 0.2 and x = 0.9 for the analytical solution. Since the CFL number is based on the 
highest velocity value (2 in this case) and is taken as 0.98, the corresponding wave 
is less diffused at a; = 0.9, contrary to the front wave located at a; = 0.2 moving at 
velocity V2 = 1. Note that in these areas, p\ = p%, since e = or e < e\. The borders of 
these areas are clearly seen at x s=s 0.28 and a; = 0.8. The constant profiles for p2 and 
pi observed respectively between a; = 0.5 and a; = 0.7 and between a; = 0.5 and a; = 0.8 
correspond to the one observed in the analytical solution. In the area between x = 0.28 
and x = 0.5, the density profiles would be expected to be constant, with the same value 
as before. Instead of that, one observes a peak value for p2 and a low value for pi, 
the sum P1+P2 being constant. This results from a coupling between the behaviour 
of the quadrature method when e tends to zero and the numerical diffusion. We are 
here in the situation q > and e small but e>e\. See Lemma 12.31 Further away from 
the discontinuities the density values tend to their analytical value. 
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Fig. 7.9. Moment dynamics for a free boundary connecting areas where e = and e>0. Results 
at i = 0.2. Top: Analytical (dashed line) and numerical (solid line) density (left) and momentum 
(right). Bottom: weights (left) and abscissas (right). The solid line corresponds to the set (p\,v\), 
and the dashed line with circles to the set (p2,V2)- 



Figure [7TU] displays the profile of quantities and at time t = 0.2. First, it 
is interesting to note that jj-^ is naturally bounded by 1 in the numerical simulation, 
as shown in Fig 17.101 such that the bound imposed with rj = 2 in never effective. Let 
us underline that in the present case, the numerical scheme allows to preserve the 
proposed cone associated with a maximal distance between the abscissas which is also 
invariant. Comparing the analytical value of |g/(Moe)| to the numerical resolution 
in Fig. 17.101 it is instructive to observe that the numerical diffusion is creating zones 
where it can be far from zero, whereas it should be zero in the analytical solution. 
However, this is done is such a way as to preserve the maximal value foreseen as the 
maximal distance between the abscissa at time t = which is one. Let us also underline 
that we have rerun this case with various values of e\ varying from 10 -9 up to 10~ 7 
without any effect on the solution. It can also be noticed that the quantity q/e 3 ^ 2 has 
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Fig. 7.10. Moment dynamics for a free boundary connecting areas where e = and e>0. 
Results at t = 0.2. Left: ratio q/(Moe). The dashed line with crosses represents the analytical 
solution, whereas the solid line represents the numerical solution, bounded by its extremal initial 
values (dotted- dashed curve). Right: ratio g/e 3 / 2 . 



large values in the regions of connection between the interior and the frontier of the 
moment space, but has no reason to be naturally limited as opposed to the previous 
quantity as it can be seen in Fig. l7.9I bottom left. 

Concerning the convergence behavior of the numerical solution, we display the er- 
ror in L\ norm relative to the analytical solution to assess convergence quantitatively, 
Tab l7.1[ for 400, 800, 1600, and 3200 cell grids. As expected, we get an experimental 
order of convergence of 0.5. These data clearly show the convergence towards the 
analytical solution for each moment. 

As a consequence, it can be seen that we have designed the proper theoretical 
setting for the transition from the interior of the moment space towards its frontier 
since the cone we have defined seems to be automatically preserved by the kinetic 
scheme we have used. Such a point would be worth a detailed study which is beyond 
the scope of the present paper. 



Grid size 


400 


800 


1600 


3200 


too 


0.049 


0.0344 


0.0244 


0.0172 


TOl 


0.0442 


0.0307 


0.0219 


0.0153 


TO2 


0.0396 


0.0271 


0.0195 


0.0136 


TO 3 


0.0362 


0.0244 


0.0177 


0.0122 



Table 7.1. L\ error on moments relative to the analytical solution. 



8. Conclusion 

In this paper, we have extended the notion of entropic measure solution of 
quadrature-based moment method for kinetic equations. Such kinetic equations are 
frequently encountered in many application fields where a complex dynamics in phase 
space is involved. Following the contribution of Bouchut & al (2003) for the pres- 
sureless gas dynamics which is the one-node quadrature version of a more general 
system of conservation laws for quadrature-based moment models, we have been able 
to provide a few problem test-cases showing that the numerical solution of the re- 
sulting system of conservation laws through kinetic schemes reproduces the defined 
entropic solution as well as the proper theoretical setting for the transition from the 
interior to the frontier of the moment space. It is an important point for the case 
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of PTC where the solution remains smooth and where the scheme allows to describe 
the phase space dynamics properly as well as for cases where the complexity of the 
dynamics in phase space leads to generalized 5-shocks, as observed for pressureles gas 
dynamics due to the weakly hyperbolic structure of the system of conservation laws. 
Two stumbling blocks still remains to be treated. First, we would need a uniqueness 
theory and a convergence analysis in a general framework in order to fully justify the 
use of the kinetic schemes for the simulation of such models. However, as explained 



order to provide uniqueness since one can exhibit multiple entropic solutions for mea- 
sure solutions. Let us underline that it is easy to construct the same type of measure 
solutions for system (|2.3p . which is the exact same collision case used by Bouchut, 
but with a motionless Dirac delta function in density localized at the collision point 
of the other two incoming "particles" . An infinite set of entropic solution can then 
be exhibited depending on the nature of the collision. As a consequence, it would 
be first useful to investigate such a point on the pressureless gas dynamics and then 
to extend it to the present system of higher order quadrature-based moment models. 
Besides, the construction of fully high order methods is still an open question and 
requires further developments. 

At last, let us mention that most of the results of the present paper do naturally 
extend to higher order moment systems, but at the price of algebra complications. In 
fact, the key point lies in the extension of the proposed study of the behavior of the 
quadrature at the frontier of the moment space, namely when the two velocities V\ 
and i>2 become equal. If we consider for instance the 6-moment model and assume 
that one of the three velocities Vi , v-i and V3 is smooth while the other two become 
equal, we are in the same framework as in the present paper. But the case when the 
three velocities tend to be equal needs to be precised in a future work. 
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